function phiN = initialShape(nelx,nely)
%Make an initial shape with several holes.
phiN = ones(nelx+1,nely+1);
r = nely/8.0;%Radius of holes
rc = [nelx/6.0 0; nelx/2.0 0; nelx*5.0/6.0 0;
     0 nely/4.0; nelx/3.0 nely/4.0; nelx*2.0/3.0 nely/4.0; nelx nely/4.0;
     nelx/6.0 nely/2.0; nelx/2.0 nely/2.0; nelx*5.0/6.0 nely/2.0];
%Center of holes
for l = 1:length(rc)
    center = rc(l,:);
    for i = 1:nelx+1
        for j = 1:nely/2+1
            if norm(center - [i j])<r
                phiN(i,j) = 0;
                phiN(i,nely-j+2) = 0;
            end
        end
    end
end